Contact information

Bruno A. Galindo

Introduction

After finish you GWAS analysis you get a list with significant associations between a marker (SNP) and your trait (phenotype), this list with from this list you will need 3 columns for the following analyzes, the SNP name, the chromosome and position of your markers.

In the image below you can find a summary of Steps, input and output information.

My Image

VARIANT EFFECT PREDICTOR (VEP)

https://www.ensembl.org/Homo_sapiens/Tools/VEP

After perform the GWAS analysis we got as output a list of significan SNPs that are associated with our trait.

The very first step to understand what this SNPs could represent can be a VEP analysis. With this analysis get more information about the location of our variants (SNPs generated by mutations) and therefor infer some consequences of them. We are gonna check if our variant falls into:

But, first we need to proivde the rsID of our variants.

Getting the rsID

The rsID is a “stable” ID that doesn’t change along the time. So is a better use than “comercial” SNP Name. We are gonna get this IDs using the website SNPchimp https://webserver.ibba.cnr.it/SNPchimp/index.php/links

Using the tab BROWSE hit the link for your specie.

My Image

Now, paste your SNPnames in the search area separated by comma.

My Image

Now you have your rsID and you can download this output in CSV or TSV format.

My Image

Running VEP

This session is still under construction!

You can find the VEP in the website https://mart.ensembl.org/info/docs/tools/vep/index.html

My Image

My Image

My Image

GALLO

Now that you have some insights about the location and consequences of your variants we can move forward and discover which genes and QTL are surrounding our significant SNPs. This step we are going to perform using GALLO package.

Many informations from GALLO you can find on (Fonseca et al. 2020)

To use GALLO you need three different input files:

Enrichment analysis

This section is based in (Dessimoz and Škunca 2017)

Now we have our genes and QTLs surrounding our significant SNPs. We can check for the presence of this genes in specific terms, and test the significance of this “presence”.

ps. from now every time I use the word term could be a GO term (GO: Molecular Function, GO: Cellular Component, GO: Biological Process), KEGG Pathway, REACTOME Pathway and other terms.

Gene enrichment using GPROFILER

Imagine we have the scenario below:

My Image

In simple terms we want to know if we have a bigger proportion of annotated genes to a specific term in our sample than in the whole population (chromosome or genome). For instance, is the proportion of genes present in the pathway grater in our sample than in the genome/chromosome? Our H_0 is what I observed in my sample, is this what we would expect just by chance?

To estimate the probability of this scenario we can use the following formula.

My Image (Source: (Dessimoz and Škunca 2017))

With this formula I’ll calculate the probability of our scenario.

But we need to check its significance and to do this we can estimate a p-value that consist in adding two probabilities under the assumption that the null hypothesis is true: 1) a scenario as extreme as our 2) scenarios more extreme than, the one actually observed.

My Image (Modified from: (Dessimoz and Škunca 2017))

In our example we found that we have 0.044 (less than 0.05 or 5%) that our scenario or more extreme scenarios happen just by chance. so we reject H0.

So we can say that our term (pathway or GO term) is enriched (more frequent / present) in our sample than in the whole genome/chromosome.

For more detailed explanation about p-value calculation you can watch the videos below:

https://youtu.be/JQc3yx0-Q9E?si=UPeDf8rAG-4u1SSA

https://youtu.be/udyAvvaMjfM?si=mY7VQLW3BbJaWFgx

Multiple testing control

As we can test for several terms (pathways) we need to control for multiple testing because, otherwise we can get many false positive just by chance.

Another way to interpret p-value is the number of false positive we “accept”, for example for a p-value of 5% we assume that we can have 5% of false positive among our significant test. Imagine a situation when you perform 400 p-values calculation, if you assume a 0,05 threshold 400 x 0,05 = 20 false positives.

To control false positive in situations of multiple testing we can use for example: Bonferroni or FDR.

In the case of FDR we will correct the p-values for the number of tests, and as consequence usualy we are going to “increase” the p-value, so after this the terms that remain positive we can consider for our next steps as significant enriched terms.

You can find extra information about FDR in this youtube video https://youtu.be/K8LQSvtjcEo?si=rEc6bRkOt0luJhiW

QTL Enrichment using GALLO

The theory behind the QTL enrichment analysis is pretty much the same that the Gene Enrichment, but instead of testing for presence in a GO term we are going to test for a specific trait, look the scenario bellow, we tried to correlate the image with the QTL enrichment output from GALLO.

My Image

  • The red circle is trait that is being tested for QTL enrichment
  • The blue circle is the chromossome we are considering our “population”, but could also be for the whole genome and not for chromosome
  • The green arrow shows the total number of QTL in chromosome 2 surrounding our significant SNP markers, this represent our sample size or n
  • The black arrow show the total # of QTLs in chromosome 2 surrounding our significant SNP markers associated with a specific trait, i.e out of the n how many are associated to a specific trait and it is the nt value.
  • The purple arrow show the total number of QTL in our population (chromossome or genome), this represents the m value
  • The red arrow show the total of QTLs in chromosome 2 associated with a specific trait, i.e out of the m how many are associated to a specific trait and it is the mt value.

The probability formula is:

\[ P = \frac{{\binom{mt}{nt} \cdot \binom{m-mt}{n-nt}}}{{\binom{m}{n}}} \]

Using our example:

\[ P = \frac{{\binom{27}{2} \cdot \binom{2967}{350}}}{{\binom{2994}{352}}} \]

#R Example
chr_ann <- 27
sample_ann <- 2
chr_qtl <- 2994
sample_qtl <- 352
log_numerator <- lchoose(chr_ann, sample_ann) + lchoose((chr_qtl - chr_ann), (sample_qtl - sample_ann))
log_denominator <- lchoose(chr_qtl, sample_qtl)
log_P <- log_numerator - log_denominator
P <- exp(log_P)

But to estimate the p-value you need to sum up all probabilities equal or rarer, in other words I need to calculate the probabilities in the same way that in the previour formula for all situations bellow and then sum up:

2 QTLs out of 27 QTLs (my sample)
3 QTLs out of 27 QTLs
4 QTLs out of 27 QTLs
5 QTLs out of 27 QTLs
.
.
.
27 QTLs out of 27 QTLs

\[ P = \sum_{i=2}^{27} \frac{{\binom{mt}{nt_i} \cdot \binom{m-mt}{n-nt}}}{{\binom{m}{n}}} \]

The output of the formula above will be the p-value and can also be calculated in r using the code below.

#R example
chr_ann <- 27 # These are the the QTLs for a specific trait annotated in a specific chromosome (mt)
chr_qtl <- 2994 # These are ALL QTLs (regardless the trait) annotated for a specific chromosome (Population size = m)
sample_qtl <- 352 # These are the SAMPLE QTLs (regardless the trait) annotated for your sample at a specific chromosome (Sample size = n)
log_denominator <- lchoose(chr_qtl, sample_qtl)

# Initialize a vector to store the results
probabilities <- numeric(26)  # We'll have 26 new combinations from 2 to 27

# Loop through different sample sizes from 2 to 27 --> these are the QTLs for a specific trait annotated in a specific chromosome for my SAMPLE (nt)
for (sample_ann in 2:27) {
  # Calculate the logarithm of the numerator
  log_numerator <- lchoose(chr_ann, sample_ann) + lchoose((chr_qtl - chr_ann), (sample_qtl - sample_ann))
  
  # Calculate the logarithm of P
  log_P <- log_numerator - log_denominator
  
  # Calculate P
  P <- exp(log_P)
  
  # Store the result
  probabilities[sample_ann - 1] <- P
}

# Print the probabilities
print(probabilities)

# Calculate the sum of probabilities
total_probability <- sum(probabilities)
total_probability

QTL enrichment

You need to set up different parameter to check the p-value for each QTL
  • qtl_db = the .gff file
  • qtl_file = the first output from Gallo with the surrounding QTLs
  • qtl_type = c(“QTL_type”, “trait”),
  • enrich_type = c(“genome”, “chromosome”),
  • chr.subset = NULL,
  • nThreads = NULL,
  • padj = c(“holm”, “hochberg”,“hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”)) 
#GALLO
#QTL enrichment analysis 
out.enrich<-qtl_enrich(qtl_db= qtl_dbb, 
                       qtl_file= out.qtl, qtl_type = "Name",
                       enrich_type = "chromosome", chr.subset = NULL, 
                       padj = "fdr",nThreads = 2)

The output has this columns

  • QTL: The QTL class or trait used for the enrichment;
  • CHR: The chromosome for that specific QTL or trait (if the option “chromosome” is informed to the argument enrich_type);
  • N_QTLs: Number of observed QTLs or traits in the dataset;
  • N_QTLs_db: Number of each annotated QTL in the qTL database;
  • Total_annotated_QTLs: Total number of annotated QTLs;
  • Total_QTLs_db: Total number of QTLs in the QTL database;
  • pvalue: P-value for the enrichment analysis;
  • adj.pval: The adjusted p-value based on the multiple test correction selected by the user;
  • QTL_type= The QTL type for each annotated trait.

QTLenrich_plot

You need to set up the following parameters
  • qtl_enrich: The output from qtl_enrich function
  • x: ID column to be used from the qtl_enrich output

  • pval: P-value to be used in the plot. If “p_value” informed, a non-adjusted pvalue will be plotted. If “p.adj” informed, the adjusted p-value from the qtl enrichment analysis will be plotted.

    Before plot the enrichment results, a new ID column will be created in order to make easier to identify the enrichment results per chromosome. Additonally, we are going to match the QTL classes for each trait and filter the top 10 enriched QTLs.

    Creating a new ID composed by the trait and the chromosome

    #GALLO
out.enrich$ID<-paste(out.enrich$QTL," - ","CHR",out.enrich$CHR,sep="")

    Match the QTL classes and filtering the Reproduction related QTLs

    #GALLO
out.enrich.filtered<-out.enrich[which(out.enrich$adj.pval<0.05),]

    Plotting the enrichment results for the QTL enrichment analysis

    #GALLO
QTLenrich_plot(out.enrich.filtered, x="ID", pval="adj.pval")

Gene Enrichment

GPROFILER

To perform the Gene enrichment analysis we are goint to use the application GPROFILER, you can use GPROFILER either in R or online.

### enriquecimento genico
#install.packages("gprofiler2")
library(gprofiler2)

#Para conferir a lista de organism -> https://biit.cs.ut.ee/gprofiler/page/organism-list

#Obs: eu entro com os ids ENSOAR...
query <- read.table ("/out_genes_60k.txt", header = T)
query <- query[,c("gene_id")]

gene_enrich <- gost(
  query,
  organism = "chircus", #o de cabra é "chircus"
  ordered_query = FALSE,
  multi_query = FALSE,
  significant = TRUE,
  exclude_iea = FALSE,
  measure_underrepresentation = FALSE,
  evcodes = TRUE,
  user_threshold = 0.10,
  correction_method = c("fdr"),
  domain_scope = c("annotated"),
  numeric_ns = "",
  sources = NULL,
  as_short_link = FALSE,
  highlight = FALSE
)


#str(gene_enrich) # para ver o formato dos meus dados

#selecionando apenas as informações da lista que me interessam para fazer meu data.frame 
result_enrich <- data.frame(gene_enrich$result)
result_enrich <- data.frame(Category = result_enrich$source,
                            ID = result_enrich$term_id,
                            Term = result_enrich$term_name,
                            adj_pvalue = result_enrich$p_value,
                            id_ensembl = result_enrich$intersection)

write.table(result_enrich,"/home/bambrozi/Sec_ARS1_GrassHill_1/GPROFILER/gene_enrich.txt", col.names=TRUE, row.names=FALSE, sep="\t", quote=F)



gostplot(
  gene_enrich,
  capped = TRUE,
  interactive = T,
  pal = c(`GO:MF` = "#dc3912", `GO:BP` = "#ff9900", `GO:CC` = "#109618", KEGG =
            "#dd4477", REAC = "#3366cc", WP = "#0099c6", TF = "#5574a6", MIRNA = "#22aa99", HPA =
            "#6633cc", CORUM = "#66aa00", HP = "#990099")
)

#to use this function to save, it is necesssary choose interactive =F in the previous step
#to use this function create a object "publish_gostplot" in the previous step.
publish_gostplot(
  gostplot_obj,
  highlight_terms = NULL,
  filename = NULL,
  width = NA,
  height = NA
)

#CREATE A TABLE

publish_gosttable(
  gene_enrich,
  highlight_terms = NULL,
  use_colors = TRUE,
  show_columns = c("source", "term_name", "term_size", "intersection_size"),
  filename = "result_table.pdf",
  ggplot = TRUE
)

GENE NETWORKS

STRING

https://string-db.org/ This session is still under construction!

Reference list

Dessimoz, Christophe, and Nives Škunca, eds. 2017. The Gene Ontology Handbook. Vol. 1446. Methods in Molecular Biology. New York, NY: Springer New York. https://doi.org/10.1007/978-1-4939-3743-1.
Fonseca, Pablo A S, Aroa Suárez-Vega, Gabriele Marras, and Angela Cánovas. 2020. “RPubs - GALLO_vignette.” https://rpubs.com/pablo_bio/GALLO_vignette.